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Abstract 

The structure of the Cu(llO) surface is studied at high temperatures using a 
combination of lattice-gas Monte Carlo and molecular dynamics methods with 

X' 

' identical many-atom interactions derived from the effective medium theory. 

The anisotropic six-vertex model is used in the interpretation of the lattice- 
gas results. We find a clear roughening transition around Tr = 1000 K, 
and Tr/Tm = 0.81. Molecular dynamics reveals the clustering of surface 
defects as the atomistic mechanism of the transition and allows us to estimate 
characteristic time scales. For the system of size 50x50, the time scale of the 
local roughening at 1150 K of an initially smooth surface is of the order of 
100 ps. 
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I. INTRODUCTION 



Thermal disorder of fcc(llO) metal surfaces has attracted considerable theoretical and 
experimental interest during the past few years. Due to its high structural and energetic 
anisotropy, this surface provides an intriguing case where a number of competing disordering 
mechanisms can be found.0'i Roughly speaking, the (110) surfaces can be divided into two 
categories. Metals such as Cu and Ni preserve their bulk-terminated (lxl) structure as the 
ground state, consisting of nearest neighbour chains separated by the conventional fee lattice 
constant. On the other hand, 5d metals Ir, Pt and Au have the tendency of forming the 
so-called (1x2) missing-row reconstructed structure or even further a (lxn) micro-faceted 
structure, where only every n:th chain is left in the first crystal layer, and the surface profile 
has (111) facets in the atomic scale. For the non-reconstructed case, we can choose copper 
as the prototype. 

Anomalous increase of disorder has been observed at Cu(110) between 500 K and 1000 
K by using X-rayi and helium diffraction! methods, low-energy ion scattering and inverse 
photoemission (IPE),i low-energy electron diffraction (LEED),i high- resolution electron- 
energy-loss spectroscopy (HREELS)i and, most recently, the impact-collision ion-scattering 
spectroscopy (ICISS).! It is now generally accepted that strongly anharmonic vibrations of 
surface atoms (anomalous Debye- Waller effect) are responsible for the increase of disorder at 
the surface above 550 K. It has been conjectured that these vibrations induce a roughening 
transition where the lattice structure of the surface is preserved, while the height-height 
correlations between the surface atoms diverge logarithmically. This was the initial conclu- 
sion from the X-ray diffraction experiments by Mochrie,! who estimated the lower bound of 
the transition temperature Tr to be 870 K. This result was later criticized by Zeppenfeld 
et alE who showed that the line shape of the specularly reflected thermal helium beam is 
not consistent with a rough surface with logarithmic height-height correlations up to 900 K. 
By extending the measurements recently up to 1200 K, they find in fact that Tr = 1070 
K.0 Also Diirr et a/.i have recently observed structural changes above 1000 K, the nature of 
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which remains unclear. It has been known for a long time that the (110) surface is missing 
from the equilibrium crystal shape of copper just below the melting temperature Tjv/,0 but 
the nature of the disordering mechanism has not been known for the temperature range 
of 1000 K up to Tm- Specifically, the possible existence of both the roughening and the 
surface premelting transitions,0 as well as their possible mutual interactions, have not been 
resolved. 

The fcc(llO) face is the most open one among the low-index faces of fee lattice, with a 
structure resembling the (more open) stepped (113), (115),... surfaces. The basic excitation 
of a stepped surface is a kink on the step edge, whereas at a low-index face, roughening is 
induced by formation of vacancies and adatoms or excitations with a still higher energy.ll! 
In the limit of very large anisotropy, the anisotropic BCSOS (body-centered solid-on-solid) 
model of the fcc(llO) face and the terrace-step-kink model of high-index surfaces become 
equivalent.0 In the case of copper (110) surface, anisotropy is not very large, and its rough- 
ening is best described by that of a low-index face. This kind of roughening has been 
described by various Monte Carlo simulations on lattice-gas and solid-on-solid modelsiil 
Recently, molecular dynamics simulations using many-atom interactions have shown pro- 
nounced disorder and finally premelting of the (110) surface of Al, Ni, Au and Cu-EHEil 
Identification of the rough phase is difficult in the molecular dynamics simulations, how- 
ever, because tractable sample sizes and simulation times preclude the identification of a 
logarithmic behaviour of the height-height correlation function. 

In this paper we use a new theoretical approach to investigate the roughening of the 
Cu(110) surface. By extracting the interatomic interactions from the well-tested effective- 
medium theory,0 the structure of the Cu(110) surface at high temperatures is studied by 
means of both the lattice-gas Monte Carlo (LGMC) and molecular dynamics (MD) simu- 
lations with comparable system sizes.S A brief account of the results was already given in 
Ref. [19], but here we shall give a detailed description of the method and of the various results 
it can provide. 

In the LGMC simulations we identify the rough phase from the behaviour of the height- 



height correlation function 

G(r) = ((h(6)-h(f)r) , (1.1) 

where h(r) is the height of the surface at point f=(x, y). Below the roughening temperature 
Tr, the interface width is finite with a finite correlation length, while at and above Tr the 
correlation function diverges logarithmically, 

G(f) ~ 2A(T)\nr +c ; T>T R . (1.2) 

At the transition point the coefficient A has the universal value A(Tr) = 1/tt 2 $ Using the 
LGMC model we find a clear roughening transition around 1000 K. The dynamics of the 
rough surface is studied via MD simulations employing the same interactions as used in the 
LGMC method. The LGMC results are qualitatively similar to the MD results which include 
the effects of lattice vibrations and relaxation. Furthermore, the atomistic mechanism of 
the roughening is found through MD to be the creation of defect clusters of vacancy and 
adatom type, and the relevant time scales can also be estimated. It is noteworthy that the 
estimated Tr is, for the first time, in a realistic temperature range in comparison with the 
previously calculated (by MD) surface premelting and bulk melting temperatures for the 
same interatomic potential.^ 

The paper is organized as follows. Section || describes the models which have been em- 
ployed, and gives a technical description of our simulations, including also a discussion of the 
effects of the inherent differences between our LGMC and MD models. The determination 



of the roughening temperature by LGMC simulations is described in Section III. In Section 



rV| we discuss the results of the MD simulations, found by following the evolution of the 
surface from a smooth to a locally rough phase. We follow in particular the time evolution 
of the occupation numbers of the surface layers and of the cluster distributions, and atomic 
diffusion. Our conclusions are given in Section |V[ 
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II. MODELS AND METHODS 



The effective medium theory (EMT) is an approximative method for calculating the 
total energy of an arbitrary arrangement of metal atoms.0 EMT relates the energy of an 
atom in metal to the local electron density and, therefore, implicitely includes many-atom 
effects. It gives a reasonable description of a number of surface properties including surface 
relaxation, reconstruction and phonon spectra. We have extended the interaction range 



in EMT to the third nearest neighbours in the fee copper lattice as outlined in Ref. [16 . 
These interactions have been used in both the lattice-gas Monte Carlo and the molecular 
dynamics simulations. 



A. Lattice-gas Monte Carlo simulations 

In our Monte Carlo simulations we use a face-centered cubic lattice-gas model with many 
body interactions derived from EMT in the form of a function E=E{Ci 1 C2 1 C^) 1 where Cj 
is the number of the j neighbours for an atom in the lattice. In Fig. [I] we show E as a 
function of C\ and C<i with C3 = 2C\. The broken curve is the derivative dE/dCi or the 
energy cost of a broken nearest neighbour bond. We also show in this figure the simplest 
possible pair interaction model (dotted line) that is normalized such as to produce the right 
cohesive energy in a perfect lattice and zero energy at zero coordination. The basic excitation 
energies at Cu(llO) are about 2.5 times higher in the pair interaction model than in the 
EMT. 

The EMT lattice-gas Hamiltonian can be written in the form 

H - fiN = Y,ME(C lu C 2i ,C 3i ) - Vi - n] , (2.1) 

where n, = 0, 1 is the occupation number of the lattice site i, Cji is the number of its 
jth neighbours anc l ^ j s the chemical potential. The external potential v is needed for the 
formation of an interface in the lattice-gas model.il We have assumed the form v, t = a for the 
first layer, and = bz~ 3 for the other layers, where is the distance of the lattice site i from 
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the 'substrate' or, equivalently, the layer number (zi = 1, 2, ..). With our choice a = 6000 K 
and b = 9000 K the formation of the first 'adsorption layer' is essentially two-dimensional 
and the bulk limit is reached for Z{ ~10. 

For comparison with our lattice gas results, we have also mapped the lattice gas interface 
to an anisotropic six-vertex solid-on-solid model or the BCSOS model. A schematic 
picture of the fcc(110) face seen from above, and its mapping to the six-vertex model, 
are shown in Fig. 0. In Fig. 0a, the black dots denote surface atoms at level and the 
open symbols with + and — are surface atoms at levels +1 and —1, respectively. The 
corresponding vertex configuration is also shown. The energy parameters of this model are 
determined from the surface excitation energies of EMT. The energy cost of broken bonds, 
ei = £2 in the nearest neighbour direction and e% = in the next nearest neighbour direction, 
are found from the energy function E = ^(Ci, C*2, C3) for a coordination corresponding to 
the uppermost atom layer for which C\ ~ 7 and C2 — 4. These energies are then scaled such 
as to reproduce the right energy for an adatom-vacancy pair on an otherwise flat face. In 
this way we find for the vertex energies: ei = e<i = 1485 K (the derivative of E with respect 
to Ci), e 3 = e 4 = 322.5 K (the distance between the curves for different C 2 in Fig. [I|), and 
65 = €q = 0. The model is in fact a pair potential model with an energy scale derived from 
a many-atom model. Notice that no external potential is needed here to define the location 
of the interface. 

The exact transition temperature Tr in the thermodynamical limit of the six-vertex 
model is given byi! 

A(T R ) = -1; A(T) = (a 2 + b 2 - c 2 )/2ab , (2.2) 

where a = exp(— ei/kT), b = exp(— e 3 /kT) and c = exp(— e 5 /kT). Using the energy pa- 
rameters given above, the roughening temperature of the EMT six-vertex model in the 
thermodynamical limit is T^(L— hdo) = 1090 K. The behaviour of the coefficient A defined 
in Eq. ( |1.2j ) is also known exactlyil as a function of temperature, 

A(T) = [ivr 2 - 71 arcsin A(T)]" 1 (2.3) 
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= 1/tt 2 + {T- T R y/ 2 [ ai + a 2 (T - T R ) + ...] . 

In Eq. (|2.3|) , the value of A at the transition point, A(T R ) = 1/ti 2 , is universal. We will 
use the leading square root (Kosterlitz-Thouless) behaviour near T R of coefficient A in 
interpreting the results of our many-body-lattice-gas Monte Carlo calculations. 

It is evident from Eq. ( |2.2| ) that the critical temperature of the six-vertex model depends 
on two parameters, the energy scale parameter e = £1 — 65, and the anisotropy of vertex 
energies v = (e 3 — e^)/{e\ — £5). In our parametrization the anisotropy ratio is v = 0.217. 



In Ref. 27 the vertex energies were determined from surface energies crxio, ctiooj &iu, and 



the anisotropy ratio was found to be v = (er m y 3/2 — 0110)/ (criooV^ — <7no). In our EMT 
lattice-gas model this method of defining the vertex energies would give v = 0.141. This 



value is still much higher than the value v = 0.080 obtained in Ref. |27] from the relaxed 
surface energies of the embedded atom method. In that work the low v value leads to a 
much lower estimate for the critical temperature. This discrepancy is mainly caused by 
the difficulty in including in the energy scales of a pair interaction model the many-body 
effects of the true interactions. Simulations with realistic many-atom interactions are thus 
really needed to settle this question. We shall show below that, just above T R , typical 
excitations at the Cu(110) surface are not large (100) or (111) facets, but monoatomic steps 
with a temperature dependent kink density. This behaviour indicates that the effects on the 
surface excitation energies of both many-atom forces and relaxation are very complicated. 

We use the method of equilibrium Monte Carlo simulations, which is knownS to be 
suitable for lattice-gas and solid-on-solid models. The many-body nature of the Hamiltonian 
Eq. (|2.1|) , and the extension of the interactions up to third neighbours in the lattice-gas 
model, increases the computation time by approximately a factor of 20 in comparison with 
the pair potential model we have also used. The lengths of the runs varied from a few 
thousand (for isotherms) to 500000 Monte Carlo steps per lattice site, for systems of size L 2 
with 12<L<196. 
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B. Molecular dynamics simulations 



The lattice-gas model discussed above lacks two major physical features. First, there is no 
real physical time connected with the evolution in the lattice. Second, the effects of lattice 
vibrations and surface relaxation are ignored. To study these effects we have performed 
large-scale molecular dynamics simulations employing the same interaction potential from 
which the site-site interaction energies are calculated for the LGMC model. This potential 
is fully documented in a previous work on surface premelting.0 

Our sample consists of 8 dynamical (110) layers on top of 4 static substrate layers, with 
either 30 close-packed [110] rows of 30 atoms or 50 rows of 50 atoms. The number of 
dynamical atoms is then 7200 and 20000 in the smaller and the larger sample, respectively. 
The dynamics is realized by the Nose-Hoover canonical thermostat.^! This method is a 
natural choice in our simulations since, throughout this work, we wish to compare the 
molecular dynamics results with those of the LGMC model, the lattice configurations of 
which are generated at constant temperature. The lattice parameter of the substrate is 
adjusted such as to correspond to the simulated temperature as defined from the thermal 
expansion produced by the interaction potential. The equations of motion are integrated 
by using the velo city- Ver let algorithm,^ which has been modified to handle the thermostat 
equations. Due to the good stability of the algorithm, we have been able to use a large 
time step of 14 fs, which is roughly a tenth of the Debye period for copper. The large 
number of atoms in MD requires special methods for selecting the pairs of atoms for which 
the interactions are calculated. The method we have used is described in the Appendix. 

The results which we shall report in this paper are those for temperatures T = 1100 K 
and T = 1150 K. In the T = 1150 K case we shall analyse a pair of runs, one initiated from 
a rough lattice-gas configuration generated by the LGMC model with a number of atoms 
corresponding to full layers (run A, 30x30 sample, length 406 ps), and the other initiated 
from an undefected smooth surface (run B, 50x50 sample, length 616 ps). At T = 1100 K 
we have made only a single run (run C, 50x50 sample, length 560 ps), initiated from the 
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undefected smooth surface. The reason for studying these temperatures is that, according 
to the LGMC results, they are within the range where the surface is rough but not yet 
premelted (the surface premelting of Cu(llO) occurs in our MD model at T > 1200 K@). 
We have also made some simulations at lower temperatures, but due to the slowing down of 
diffusion, the statistics for the defect clusters becomes poor. 



C. Inherent differences of the LGMC and MD models 

In our LGMC simulations the energy function E is derived from the effective medium 
theory by using the bulk lattice constant at zero temperature. In the lattice-gas model there 
is no surface relaxation, but its effect on the transition temperature can easily be estimated 
by comparing the energy of the basic excitation of the unrelaxed lattice-gas model with that 
at the relaxed surface. The inward relaxationlll means higher electron density or effectively 
higher coordination, and thus decrease in the slope of the energy function of Fig. [[]. In our 
lattice-gas model we therefore overestimate the excitation energies and, consequently, the 
roughening temperature of Cu(110), by about 5%. This estimate will be used in determining 
our final estimate for Tr. 

In the LGMC energy function the effect of thermal expansion is not included, but this 
effect is in most part included in our estimate of the relaxation effects on the roughening 
temperature given above. A more fundamental difference between the two models is the 
fluctuating particle number in the LGMC simulations vs. the fixed particle number in the 
MD simulations. In practice this difference leads to a much slower sampling of uncorrelated 
lattice configurations in the MD simulations. In the MD simulations the interaction potential 
has been smoothened at the cutoff radius (between the 3 rd and 4 th neighbours) as described 



in the Appendix of Ref. [U| This has been done in such a way that effectively (through 
indirect interactions) increases the anisotropy, and thus lowers the roughening temperature 
of our MD model. We only note here that the way of cutting the interaction range in MD 
simulations can have a considerable effect on such surface properties that are very sensitive 
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to the strength of the interaction beyond the nearest neighbours. This is evident from the 
discussion of the missing-row reconstruction given in Ref. 121 . 



III. MC RESULTS 

Calculated height-height correlation functions for the lattice-gas system with L = 96 and 
for the six-vertex system with L = 50 are shown in Fig. Indicated symbols denote the 
Monte Carlo results and solid lines are best logarithmic fits to the data points. The distance 
r is in the direction of the next-nearest neighbours in the surface plane. In the lattice-gas 
model the value of the chemical potential was chosen such as to produce an interface at a 
level Zi > 10, i.e. Vi < 10 K. For temperatures above T = 1000 K, the effect on the surface 
structure of the substrate potential difference between adjacent layers is negligible. As noted 



in Ref. |3T, finite size effects become significant for r ~ L/5. On the other hand, a few lowest 
values of r have been discarded in the fits because Eq. ( |1.2| ) becomes exact only at large r. 
We have also done simulations for bigger systems, but increasing computational cost makes 
it impossible to sample the long-wavelength fluctuations properly. This finite time effect 
prevents us from being able to make fits to the correlation functions at larger values of r. 

From the logarithmic fits to the correlation functions we can infer the value of the 
coefficient A defined in Eq. ( |1.2|) . The result is shown in Figs, ^a and |}d as a function 
of temperature for the EMT lattice-gas model and for the six-vertex model, respectively. 
A square root fit of Eq. ( |2.3|) to the Monte Carlo data is shown by the solid line. From 
these fits we find T 1 ^ = 1000 K and = 1030 K. It is evident that the dependence on the 
system size of the transition temperature is beyond our accuracy. The difference between 
the exact value in the thermodynamic limit of the six-vertex model (Tjj (L— kxs) = 1090 K) 
and our numerical result (Tr = 1030 K) is caused by two effects, the effect of finite size and 
(mainly) the relatively short distances over which the fitting of the correlation function has 
been made. A systematic deviation of similar magnitude (~ 60 K) is expected to be present 
also in our lattice-gas result, where MC runs of comparable length have been performed, 
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and logarithmic fits to the correlation functions have been made in a similar way. Thus we 
can conclude that T%(L-xx) w 1060 K. 

For comparison we have also determined the Tr of our lattice-gas model by studying the 
finite size behaviour of the interfacial width , 5h 2 = ((h—h(r)) 2 )^ t , where h is the nominal 
surface height for a given configuration.ilil Above the transition temperature, Sh 2 diverges 
logarithmically as a function of the linear size of the system size, L, or 6h 2 ~ A(T) InL for 
large L. In Fig. |]c we show the result for the transition temperature as determined from 
the coefficient A(T) for the interfacial width. Here data for L = 12, 18, 24 and 30 has been 
used. The result is found to be consistent with that obtained from the correlation functions. 

The Hamiltonian Eq. ( |2.1|) , leads us to still another way of deducing the roughening 
transition, namely through layering transitions in an adsorption system.0 For a thickening 
adsorption film the critical points of the layering transitions approach the bulk roughening 
temperature.il In Fig. [| we show a few adsorption isotherms corresponding to the formation 
of the first 'adsorption layer' and the inverse slope of the adsorption isotherms at the layering 
transition point, fi — //i, as a function of temperature. The solid line in Fig. [5]b is a fit with 
the two-dimensional critical exponent of the susceptibility, 7 = 7/4, and with a critical 
temperature 7\ = 845 K. A similar analysis for the next layer gives T 2 = 870 K. We do not 
intend to determine Tr by repeating this procedure for further layers; we only note here that 
these critical layering points, T^, are lower bounds for the roughening temperature.il In our 
lattice-gas model the first adlayer grows essentially two-dimensionally (the second layer is 
almost unoccupied), which means that 7\ pa T^ d and we have T l 2 9 d /T l i = 0.80. By using vertex 
energies e± and 63 as the energy parameters in the two-dimensional Ising modelj^l we obtain 
T ising/ T 6v = q.80. Besides being a further confirmation of our estimate of the roughening 
temperature in the lattice-gas model, this result shows that a simple pair-potential model 
with relatively short-range interactions (such as the six- vertex model) can be used to describe 
the roughening transition of a metallic system, provided the values of the energy parameters 
of the model are determined in an appropriate way. The effect on the roughening transition 
of the broken particle-hole (adatom-vacancy) symmetryillzl is too small to be detected in 
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our lattice-gas model: the energy difference between convex and concave corners on Cu(llO) 
is less than 10 K. 

In Fig. H snapshots of typical configurations of a lattice-gas system with L = 30 are 
shown for two temperatures. Atoms in the lowermost odd layers are coloured black to guide 
the eye. Below T = 1000 K excitations are mainly single adatoms and vacancies at the 
surface. Around T — T R the connectivity of the adatom and vacancy clusters increases, and 
at higher temperatures the surface is clearly 'rough'. The anisotropy of the clusters reflects 
the high anisotropy of the excitation energies in the surface plane, but there is no clear 
indication of (lll)-type facetting since the anisotropy in our model is not very large. In the 
snapshot at T = 1100 K two overhangs can be seen: in both cases one of the four nearest 
neighbours at the next atom layer below is missing. Near T = T R the overhang density is 
very low (<0.001), and the solid-on-solid approximation of the six-vertex model is valid. 

We conclude that the roughening temperature of our EMT lattice-gas model is T l R 9 = 
1060 K. Considering the effect on the energy scale of surface excitations of relaxation (see 
Sec. [11 C| ), our best estimate for the roughening temperature of Cu(110) is T R « 1000 K. 
Comparing this with the bulk melting point for the same potential, Tm = 1240 K,@ we find 
T R /T M = 0.81. 

IV. MD RESULTS 

We shall first compare the runs A (rough initial configuration) and B (smooth initial 
configuration) by studying the occupation numbers of the surface layers as a function of time. 
Atoms are attached to layers according to their z-coordinate by dividing the z-dimension of 
the sample in slices of thickness Az = a/2y / 2, where a is the conventional fee cubic lattice 
constant. Az then corresponds to the distance between adjacent (110) planes in the bulk. 
We adopt hereafter the convention of 'crystal' and 'adatom' layers. The crystal layer c\ is 
the uppermost full layer of the undetected surface at T = 0, and the adatom layer a\ is the 
first unoccupied layer at T = 0. The other crystal and adatom layers are numbered from 
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the T = interface into the bulk and into the vacuum, respectively. 

The occupation numbers of the surface layers in runs A and B for T = 1150 K are 
plotted in Fig. |7| as a function of time. We can extract from Fig. |7] several important 
features of surface roughness. First, the initial roughness (i.e. occupation numbers) of the 
LGMC surface in the run A remains essentially the same also during the MD phase. The 
width of the interface seems to slightly increase after the MD simulation is started, which 
can be seen from the systematic increase in the occupation of the two adatom layers. The 
rough initial LGMC configuration is thus stable against the lattice vibrations and relaxation 
effects, i.e., it has no tendency of smoothing during the MD simulations. Second, we can 
reach the roughness of the run A during the first 150 ps of run B, which is initiated from 
the smooth surface. Notice that this is directly an estimate for the physical time needed 
for the surface to 'roughen' locally via its intrinsic dynamics. Third, one can see irregular 
fluctuations of a few tens of ps, on top of the normal statistical fluctuations. These long- 
wavelength fluctuations can especially be seen in the run A, which is made for the smaller, 
30x30, sample. These fluctuations are connected with the diffusion dynamics of adatom 
and vacancy clusters. They show that, during our simulation time, it is possible to produce 
at least a few uncorrelated lattice configurations, i.e., we are able to travel through the 
configuration space of the rough surface. 

At T = 1100 K (run C) the behaviour of the occupation numbers is qualitatively similar 
to that shown in Fig. [7[ However, it takes much more time (about 300 ps) to reach the 
average LGMC level. We can interpret this as a 'critical slowing down' effect, which makes 
the observation times for the cluster dynamics in MD simulations extremely long when the 
transition temperature is approached. Some snapshots of two-dimensional adatom clusters in 
layer a,\ at 1150 K (run B) are shown in Fig. |8|. The configurations in the first two snapshots 
(with a time interval of 14 ps or 1000 MD time steps) are clearly strongly correlated. In the 
third snapshot we obviously see a new lattice configuration. Despite the large vibrations 
and atomic diffusion, the clusters have a well-defined (110) surface symmetry. 

In order to study atomic mobility we have calculated the probability distributions 
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p(x, y, z; t) for the surface layers, where p is the probability that an atom is in the posi- 
tion (x,y,z) at time t, when it is at the origin at t — 0. The function p(z;t) for layer a\ 
shown in Fig. ^| illustrates the large number of layer changes during the simulations: After 
the observation time of 70 ps it is about equally probable to find the atom in layer c\ as in 
layer a%. This result is be related to the initial rising time of the occupation of the layer 
a\ shown in Fig. [7[ The p(x; t) distribution for the cross-channel [001] direction is shown 



in Fig. [1C]. The spacing between the maxima in p(x; t) corresponds to a half of the lattice 
constant, which indicates that the cross-channel [001] diffusion is the main mechanism of 
layer changes. This result has also been found from previous MD simulations of the fcc(110) 
surface diffusion.!! 

We have studied the residence time of an atom at the lattice sites of each surface layer by 
analyzing the atomic displacements during the last 42 ps of the run B (T = 1150 K), and by 
collecting the statistics of jumps from one lattice site to another. The collected distributions 



of the residence times for the surface layers are shown in a semilogarithmic scale in Fig. [TT 
After an initial rise all distributions seem to obey an exponential law, p oc e - */ 1 ", from which 
the mean residence times, r, for the diffusing atoms can be estimated to be 1.5 ps , 2.2 ps, 
3.5, ps and 4.5 ps for 02, a\, ci, and C2, respectively. These residence times correspond to 
10-30 lattice vibrations. Notice, however, that in layer C2 27% of the atoms are not diffusing 
at all during our observation time. The 'true' residence time for that layer must therefore 
be longer than the 4.5 ps quoted above. The corresponding fraction of immobile atoms in 
C\ is only 2%, and practically all atoms in layers a<i and a\ have made at least one jump 
during our observation time. 

The diffusion constants D x ,D y for the surface layers are determined from the time- 
dependent mean-squared atomic displacements (msd) r 2 (t) = ((r(£) — r(t )) 2 ), where the 
brackets mean averages taken over 20 initial times t and over atoms belonging initially to 
a given layer. D V1 v = x,y, can be obtained from the Einstein relation as, for t — > 00, 
r v = 2LU+const. We find that (D x , D y ) for a 2 ,...,c 2 are (1.80, 4.13), (2.05, 2.69), (1.71, 
2.77), and (0.88, 1.24) xl0~ 5 cm 2 /s, respectively. It is evident that the calculated surface 
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mobility is high, near the mobility level in the bulk liquid at the melting point .0 However, 
the persistence of the lattice structure makes the diffusion clearly anisotropic. We emphasize 
that what has been considered here is the intrinsic (tracer) diffusion constant, not the mass 
transport constant which is expected to become isotropic at T — Tr due to the macroscopic 
step diffusion.0il 

The time scale of the roughening of a smooth surface can also be measured by studying 



the time evolution of cluster size distributions. In Fig. |12| a we show the size distribution 
at T = 1100 K of adatom clusters in layer ai as a function of time (MD run C). After a 
simulation of 0.56 ns we are still far from the equilibrium distribution for clusters of size N > 
60. At T = 1150 K the size distribution is stable (within our statistical accuracy) after 0.16 
ns for N up to iV ~ 120. The result for the time interval 160 ps to 616 ps (run B) is shown 



in Fig. |r2|b. The dotted line is the behaviour of the lattice-gas system with of same size and 
at the same temperature. The results for the vacancy clusters in layer c\ are similar. The 
height-height correlation function G (not shown here) of Eq. ( |1 . 1| ) as determined from the 
MD data for T = 1150 K, is an increasing function of r only for r < 5. With r taken to be 
in the direction of the next-nearest neighbours (the direction of weaker bonds), this result 
reflects the fact that, during this simulation of length 0.6 ns and well above the roughening 
temperature, we have been able to collect reasonable statistics for defect clusters only for 
iV < 100. The dynamics of the cluster size distributions in our MD simulations is consistent 



with the results for the kinetics of the roughening of a stepped surface discussed in Ref. £40]: 
roughness sets in over short distances and then spreads asymptotically as t 1//3 . The time 
scale of roughening, 100 ps (Fig. |7]) at T — 1150 K and 300 ps at T = 1100 K, defined here 
as the time required for the occupation numbers of the innermost layers to saturate to their 
'rough' values, clearly describes the local behaviour. For systems of size L 2 with L > 50, this 
quantity should depend only weakly on L. We have not tried to determine the asymptotic 
behaviour for late times. 

Based on the calculated diffusion constants in the surface region, we can conclude that 
for an atom found initially in layer c±, it takes at T = 1150 K about 30 ns (50 x our 
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simulation time) to travel through a distance similar to the length of the simulation box 
in the shorter y-direction. The corresponding time for diffusion across the box in the x- 
direction is about 90 ns. These characteristic times provide an order-of-magnitude estimate 
for the amount of computation time needed for a reliable calculation of the height-height 
correlation function at one temperature by the MD method. In order to be able to determine 
Tr, a few temperatures below 1150 K should be studied. The computation time required 
by the present MD model, with the EMT interactions for Cu(110), would be of the order of 
one year of Cray X-MP CPU time. 



V. CONCLUSIONS 

We have studied the roughening of the Cu(110) surface via atomistic simulations using 
both the Monte Carlo and molecular dynamics methods which employ the effective-medium 
theory as the interaction potential. The Monte Carlo simulations made for our lattice-gas 
model, adjusted by corrections arising from relaxation and finite size effects, show clearly 
that the Cu(110) surface has a roughening transition around T = 1000 K, about 200 K below 
the surface premelting temperature ,0 determined for the same potential. The rough phase 
has been identified from the logarithmic behaviour of the height-height correlation function 
and of the interface width. The transition temperature is obtained from the temperature- 
dependent coefficient of the logarithmic term. The finite size effects have been estimated by 
mapping the lattice-gas interface to an anisotropic six- vertex model, for which the behaviour 
in the thermodynamical limit is known exactly. Knowing the bulk melting point for the 
same potential, we find Tr/Tm = 0.81. This is in excellent agreement with the experimental 
result Tr/Tm = 0.79 given in Ref. |2|. Molecular dynamics simulations of the same surface 
cannot produce statistics enough for the height-height correlation function, but they show 
the stability of the LGMC surface against lattice vibrations and relaxation, i.e., the rough 
LGMC surface shows no tendency of smoothing when used as an initial configuration for 
MD simulations. Furthermore, MD simulations show that the roughening mechanism is 
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connected with the dynamics of diffusive adatom and vacancy clusters and gives information 
about the corresponding time scales. For the system of size 50x50 studied in this work, an 
initially smooth surface 'roughens' locally in the time scale of about 100 ps at T = 1150 
K, and of about 300 ps at T = 1100 K. The time interval between two uncorrelated rough 
configurations in MD simulations is found to be of the same order. 



ACKNOWLEDGMENTS 

We wish to thank Professor Robert Swendsen for useful suggestions. This work has been 
partly supported by the Academy of Finland (JM) and by the Emil Aaltonen Foundation 
(HH). Part of the computations have been performed on Cray X-MP at the Center of 
Scientific Computing (CSC) in Espoo, Finland. 



APPENDIX: CELL METHOD AND NEIGHBOUR LIST 

A well-known problem in MD simulations with short-range interactions is the question 
of how to select only those pairs of atoms which are within the range of the potential. The 
brute-force testing of all the ~N(N — 1) pairs at each time step is a prohibitive waste of time 
even in supercomputers, if N is of the order of few thousands atoms. In our simulations, we 
have constructed the so-called neighbour list ,0 which includes the indices j of all neighbours 
of atom i for which < r^, the list radius. By using this list, only ^NNl interactions 
have to be calculated, where Ni is the average number of neighbours for a given atom. 
The cutoff rc in the EMT interaction potential! is between the third and fourth nearest 
neighbours in the T = fee copper lattice. We have set to correspond to the distance to 
the middle of the sixth and seventh neighbours, which means that Nl falls within the range 
of 80-90. In order to minimize the frequency of the updating of the list, we study at each 
time step the displacements of atoms from the previous updating step. The list is updated 
whenever the maximum displacement exceeds half of the list 'skin', tl — rc- 
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For the sample sizes used in this work, also the N 2 computations needed for updatings 
become significant. We have managed to speed up the updating using the following method. 
We divide the sample into n x x n y subcells according to the number of the 2D lattice sites, 
i.e., 30x30 or 50x50 for the samples studied in this work. For each subcell Cj, the indices 
of neighbouring subcells Cj are stored in a map in the beginning of the simulation. After 
attaching a given atom to its subcell, it is sufficient to go through its own subcell and (in our 
case) only 15 neighbouring subcells to construct its neighbour list. This method, being linear 
with the system size, clearly overrides the N 2 searches over all pairs, and brings the time 
needed for updating comparable to that needed in the calculation of the EMT interactions. 
Our method can be regarded as a 2D modification of the so-called linked-cell method,0 
with the important difference that the size of the subcell is chosen to accommodate only 
one lattice site in two dimensions. In the present method - remember that we also use a 
stable integration algorithm which allows a large time step - we have been able to perform 
MD simulations for 30000 atoms (20000 of which dynamical) on a Decstation 5000/200 (3.7 
Mflops) scalar machine at a rate of about 20 ps/day. 
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FIGURES 

FIG. 1. Energy of an atom in EMT lattice-gas (full curves). The broken curve is the derivative 
of E with respect to Ci, and the dotted line is the energy in a simple pair potential model. 

FIG. 2. Mapping of the Cu(llO) surface to the six-vertex model, (a) Correspondence between 
the fcc(llO) face and the six- vertex model (a model configuration with a vacancy on an otherwise 
flat surface), and (b) energies of the six vertices allowed by the solid-on-solid restriction. 

FIG. 3. Height-height correlation functions for (a) lattice-gas model and (b) six-vertex model. 

FIG. 4. Coefficient A as a function of temperature (a) from the correlation functions of the 
lattice-gas model, (b) from the correlation functions of the six-vertex model and (c) from the 
interfacial width of the lattice-gas model. 

FIG. 5. (a) Formation of the first 'adsorption layer' in the lattice-gas model and (b) inverse 
slope of the adsorption isotherms at the layering transition point. 

FIG. 6. Snapshots of lattice-gas configurations for two temperatures. 

FIG. 7. The occupation of surface layers 02, C2 as a function of simulation time in MD runs 
A and B for T = 1150 K. Run A was initiated from a rough lattice configuration and run B from 
an undefected smooth surface. 

FIG. 8. Snapshots of two-dimensional clusters in layer a\ in MD for T = 1150 K. 

FIG. 9. The probability distribution p(z;t) showing the atomic mobility perpendicular to the 
(110) surface. 

FIG. 10. The probability distribution p(x;t) showing the atomic mobility in the [001] direction 
parallel to the (110) surface. 
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FIG. 11. The natural logarithm of the residence time distribution for the surface layers. Shown 
also are the linear fits, from which the mean residence times r for each layer have been estimated. 

FIG. 12. Size distribution of two-dimensional clusters in layer a\ (a) as a function of time in 
the MD run C for T = 1100 K and (b) in the MD run B for T = 1150 K. The dotted line in (b) is 
the LGMC behaviour at the same temperature. 
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